# R script for for replicating main analysis
# in 'When Parties Move to the Middle: The Role of Uncertainty'
# by J. Lindvall, D. Rueda, and H. Zhai
# this file written by: H Zhai (2022-11-03 [updated: -])
# on device: Mac Pro 13 Dual-Core Intel Core i5 2.3 GHz 

# PLEASE MAKE SURE ALL THE REPLICATION FILES (DATA AND SCRIPTS) ARE STORED AT THE SAME LEVEL IN THE SAME DIRECTORY
# OR MAKE SURE THE DIRECTORY-RELATED CODES ARE PROPERLY ADJUSTED 
# TO ENSURE THE CODES RUN WITHOUT DIRECTORY-RELATED PROBLEMS
# RESTART R SESSION BEFORE RUNNING

# This file replicates results from all main analysis reported in the main body
# In the order they appear in text

# BEGIN SCRIPT
rm(list = ls())

# pkgs --------------------------------------------------------------------

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("scales")) install.packages("scales")
if (!require("lmtest")) install.packages("lmtest")
if (!require("sandwich")) install.packages("sandwich")
if (!require("stargazer")) install.packages("stargazer")
if (!require("margins")) install.packages("margins")

# load data ---------------------------------------------------------------

# make sure the .RData file is either downloaded from Dataverse or `R_0_2_make_data.R` has been run

load("RData_data_core.RData") # use core data only
load("RData_data_core2.RData") # use core data + for party-level analysis

# set graphx --------------------------------------------------------------

theme_base <- theme_get() # save for reset at the end
theme_set(theme_minimal(base_size = 14))

# Figure 1 ----------------------------------------------------------------

## define helpers ----

### calculate median measures
#### adapted from open-source contribution: https://stackoverflow.com/questions/18887382/how-to-calculate-the-median-on-grouped-dataset
#### by user 'A5C1D2H2I1M1N2O1R2T1' (2013-09-22)
gmed <- function(x, f, out = NULL) { # function for group median measures
  low <- min(x) # lower boundary of distribution
  up <- max(x) # upper boundary of distribution
  mid <- zoo::rollmean(x, k = 2) # get interval boundaries 
  while (isTRUE(all.equal(max(mid), up))) {
    mid <- mid[1:length(mid)-1] # drop duplicate boundaries
  }
  breaks <- unique(c(low, mid, up)) # boundaries
  int <- cut(x, breaks = breaks, include.lowest = T) # cut into intervals
  int <- sapply(strsplit(gsub("\\[|\\]|\\(|\\)", "", int), ","), as.numeric) # trim & convert interval to numerical
  cf <- cumsum(f) # cumulative frequency
  midrow <- findInterval(max(cf)/2, cf) + 1 # median class
  l <- int[1, midrow] # lower class boundary of median class
  h <- diff(int[, midrow]) # width of median class
  f2 <- f[midrow] # frequency of median class
  cf2 <- cf[max(midrow - 1, 1)] # cumulative frequency up to median class
  n_2 <- max(cf)/2 # total frequency divided by 2
  # return output as required
  if (!is.null(out)) {
    if (out == "med") unname(l + (n_2 - cf2)/f2 * h) # return median
    else unname(f2 / (h*max(cf))) # return class density
  } 
} 

### add bounds 
f_addbounds <- function(data) {
  p <- with(data, c(2*p[pt=="L"] - p[pt=="C"], 2*p[pt=="R"] - p[pt=="C"]))
  add_row(data, p, v = c(0,0))
} # add upper/lower bounds 

### simulate position
#### NOTE: NOT USED IN THE LATEST VERSION
f_simp <- function(x) { 
  tx <- length(x)
  p <- c(
    p[pt=="L"] + x,
    rep(p[pt=="C"], tx),
    p[pt=="R"] - x
  )
  v <- rep(v, each = tx)
  pt <- rep(pt, each = tx)
  st <- rep(x, 3)
  df <- tibble(p, v, pt, st)
  df %>% group_by(st) %>% 
    arrange(p) %>% 
    mutate(pmed = gmed(p, v, out = "med"),
           dmed = (gmed(p, v, out = "den") %>% 
                     rescale(from=c(0,1), to=c(0,100)))
    ) %>% ungroup()
}

### simulate vote share
f_simv <- function(x) { 
  tx <- length(x)
  p <- rep(p, each = tx)
  v <- c(
    v[pt=="L"] - x/2,
    v[pt=="C"] + x,
    v[pt=="R"] - x/2
  )
  pt <- rep(pt, each = tx)
  st = rep(x, 3)
  df <- tibble(p, v, pt, st)
  df %>% group_by(st) %>% 
    mutate(pmed = gmed(p, v, out = "med"),
           dmed = (gmed(p, v, out = "den") %>% 
                     rescale(from = c(0,1), to = c(0,100)))
    ) %>% ungroup()
}

## setup ----

### 3 parties
pt <- c("L", "C", "R")

### base position = {40, 50, 60}
p <- c(40, 50, 60)

### base vote = {40, 20, 40}
v <- c(40, 20, 40)

## Left Panel: Benchmark ----

### make data
base <- tibble(pt, p, v) %>% arrange(p) %>% 
  mutate(pmed = gmed(p, v, out = "med"),
         dmed = (gmed(p, v, out = "den") %>% 
                   rescale(from=c(0,1), to=c(0,100))))

### benchmark density
dmed0 <- unique(base$dmed)

### plot result 
ggplot(data = base, aes(x = p, y = v)) +
  geom_bar(stat = "identity", width = 5, alpha = 0.7) +
  geom_smooth(data = f_addbounds(base), method = "loess", color = "black", size = 0.3) +
  geom_text(aes(x = p, y = v/2, label = pt), size = 7) +
  xlim(20, 80) + ylim(0,60) + scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(title = "Benchmark Scenario", x = "Ideology", y = "Vote Share (%)") 

## Middle and Right Panels: Bimodality & Unimodality ----

### simulate vote share change
simv <- f_simv(-15:20)

### Middle Panel
ggplot(data = simv[simv$st==-15,], aes(x = p, y = v)) + # max divergence
  geom_bar(stat = "identity", width = 5, alpha = 0.7) + 
  geom_smooth(data = f_addbounds(simv[simv$st==-15,]), method = "loess", color = "black", size = 0.3) +
  geom_bar(data = base, stat = "identity", fill = NA, color = "darkgrey", lty = "dashed", width = 5) + # benchmack case
  geom_smooth(data = f_addbounds(base), method = "loess", color = "darkgrey", lty = "dashed", size = 0.3) +
  geom_text(aes(x = p, y = v/2, label = pt), size = 7) +
  xlim(20,80) + ylim(0,60) + scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(title = "Vote Share Change I: Bimodality", x = "Ideology", y = "Vote Share (%)") 

### Right Panel
ggplot(data = simv[simv$st==20,], aes(x = p, y = v)) + # max convergence
  geom_bar(stat = "identity", width = 5, alpha = 0.7) + 
  geom_smooth(data = f_addbounds(simv[simv$st==20,]), method = "loess", color = "black", size = 0.3) +
  geom_bar(data = base, stat = "identity", fill = NA, color = "darkgrey", lty = "dashed", width = 5) + # benchmack case
  geom_smooth(data = f_addbounds(base), method = "loess", color = "darkgrey", lty = "dashed", size = 0.3) +
  geom_text(aes(x = p, y = v/2, label = pt), size = 7) +
  xlim(20,80) + ylim(0,60) + scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(title = "Vote Share Change II: Unimodality", x = "Ideology", y = "Vote Share (%)") 

# clean env.
rm(list = ls()[!grepl("(data_core|theme_base)",ls())])

# Figure 2 ----------------------------------------------------------------

## load data ----
load("RData_measure_voter.RData")

## subset data ----
IDs.ind <- c("Australia", "Austria", "Belgium", "Canada", "Denmark", "Finland", "France", "Germany", 
             "Iceland", "Italy", "Luxembourg", "Netherlands", "New Zealand", "Norway", "Portugal", 
             "Spain", "Sweden", "Switzerland", "United Kingdom", "United States")
measure_voter <- measure_voter %>% 
  filter(eyear>=1965 & countryname %in% IDs.ind)

## Panel A: Median Position ----
measure_voter %>%
  group_by(eyear) %>%
  mutate(leftright_med_mean = mean(leftright_med, na.rm = T)) %>% 
  ggplot(aes(x = eyear)) + 
  geom_smooth(aes(y = leftright_med), color = "black", fill = "grey", size = 0.5, lty = "dashed") +
  geom_line(aes(y = leftright_med_mean), color = "gray30") +
  labs(x = "", y = "Median Voter Position", title = "Median voter positions over time", subtitle = "Rich democracies, 1965-2018") +
  scale_x_continuous(breaks = scales::pretty_breaks(n=10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n=10))

## Panel B: Median Certainty ----
measure_voter %>%
  group_by(eyear) %>%
  mutate(leftright_den_rlog_mean = mean(leftright_den_rlog, na.rm = T)) %>% 
  ggplot(aes(x = eyear)) +
  geom_smooth(aes(y = leftright_den_rlog), color = "black", fill = "grey", size = 0.5, lty = "dashed") +
  geom_line(aes(y = leftright_den_rlog_mean), color = "gray30") +
  scale_x_continuous(name = "", breaks = scales::pretty_breaks(n=10)) +
  scale_y_continuous(name = "Median Voter Certainty (Log)", 
                     breaks = scales::pretty_breaks(n=10)) +
  labs(title = "Median voter certainty over time", subtitle = "Rich democracies, 1965-2018") 

# Figure 3 ----------------------------------------------------------------

## load data ----
load("RData_measure_party.RData")

## subset data ----
measure_party <- measure_party %>% 
  filter(eyear >= 1965 & countryname %in% IDs.ind)

## make plot ----
ggplot(data = measure_party) +
  geom_density(aes(leftright_ml, fill = "ml", color = "ml"), # main left measure
               alpha = 0.7) +
  geom_density(aes(leftright_mr, fill = "mr", color = "mr"), # main right measure
               alpha = 0.7) +
  scale_color_grey(start = 0.8, end = 0.2, name = "Party Type", labels = c("Main Left", "Main Right")) + 
  scale_fill_grey(start = 0.8, end = 0.2, name = "Party Type", labels = c("Main Left", "Main Right")) +
  xlim(0, 100) +  # range 0-100
  labs(x = "", y = "Density", title = "Distribution of main party left-right positions", 
       subtitle = "Rich democracies, 1965-2018") 

# Figure 4 ----------------------------------------------------------------

## reshape data (for Panels 3-4) ----
ts.party  <- data_core2 %>% 
  select(eyear, main_left, main_right) %>% 
  mutate(eyear = as.integer(as.character(eyear))) %>% 
  group_by(eyear) %>% 
  mutate(
    main_left_m = mean(main_left, na.rm=TRUE), 
    main_right_m = mean(main_right, na.rm=TRUE)
  ) %>% 
  ungroup()

## Panel A: Average Position ----
measure_party %>%
  group_by(eyear) %>% 
  mutate(leftright_mid_mean = mean(leftright_mid, na.rm = T)) %>% 
  ggplot(aes(eyear)) +
  geom_smooth(aes(y = leftright_mid), color = "black", fill = "grey", lty = "dashed", size = 0.5) +
  geom_line(aes(y = leftright_mid_mean), color = "gray30") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(x = "", y = "Average Position ( [Left + Right] / 2)", 
       title = "Main party average positions over time", 
       subtitle = "Rich democracies, 1965-2018")

## Panel B: Distance ----
measure_party %>%
  group_by(eyear) %>% 
  mutate(leftright_dif_mean = mean(leftright_dif, na.rm = T)) %>% 
  ggplot(aes(eyear)) +
  geom_smooth(aes(y = leftright_dif), color = "black", fill = "grey", lty = "dashed", size = 0.5) +
  geom_line(aes(y = leftright_dif_mean), color = "gray30") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(x = "", y = "Policy Distance (|Left - Right|)", 
       title = "Main party policy distances over time", 
       subtitle = "Rich democracies, 1965-2018")

## Panel C: Average Position - Main Left ----
ggplot(ts.party, aes(eyear)) +
  geom_smooth(aes(y = main_left), color = "black", fill = "grey", lty = "dashed", size = 0.5) +
  geom_line(aes(y = main_left_m), color = "gray30") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10), limits = c(30,66)) +
  labs(x = "", y = "Average Position ( Main Left)", 
       title = "Main left positions over time", 
       subtitle = "Rich democracies, 1965-2018")

## Panel D: Average Position - Main Right ----
ggplot(ts.party, aes(eyear)) +
  geom_smooth(aes(y = main_right), color = "black", fill = "grey", lty = "dashed", size = 0.5) +
  geom_line(aes(y = main_right_m), color = "gray30") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10), limits = c(30,66)) +
  labs(x = "", y = "Average Position (Main Right)", 
       title = "Main right positions over time", 
       subtitle = "Rich democracies, 1965-2018")

# Table 1 -----------------------------------------------------------------

## define helper (clustered HCSE) ----
chcse <- function(model) { 
  coeftest(model, vcov. = vcovHC(
    model, type = "HC0", cluster = "countryname"))[,2] # hcse clustered by country (only hc0 works)
}

## fit models ----
### cols 1-3: fixed-effects (FE) models
#### col 1 
lm_mid_fe_base <- lm(data = data_core, leftright_mid~leftright_den_rlog_lag*leftright_med_lag+countryname+eyear)
lm_mid_fe_base_se <- chcse(lm_mid_fe_base)
#### col 2
lm_mid_fe_str <- update(lm_mid_fe_base, formula. = .~.+gov_party_y3+ud_ipol_y3+openc_y3+realgdpgr_y3+unemp_y3)
lm_mid_fe_str_se <- chcse(lm_mid_fe_str)
#### col 3
lm_mid_fe <- update(lm_mid_fe_base, 
                    formula. = .~.+
                      gov_party_y3+ud_ipol_y3+openc_y3+realgdpgr_y3+unemp_y3+
                      compl_leftright+compr_leftright+vturn_lag+leftright_pl_lag+effpar_ele_lag)
lm_mid_fe_se <- chcse(lm_mid_fe)
### cols 4-6: lagged-dependent-variables (LDV) models
#### col 4
lm_mid_ldv_base <- lm(data = data_core, leftright_mid~leftright_den_rlog_lag*leftright_med_lag+leftright_mid_lag+eyear)
lm_mid_ldv_base_se <- chcse(lm_mid_ldv_base)
#### col 5
lm_mid_ldv_str <- update(lm_mid_fe_str, formula. = .~.-countryname+leftright_mid_lag)
lm_mid_ldv_str_se <- chcse(lm_mid_ldv_str)
#### col 6
lm_mid_ldv <- update(lm_mid_fe, formula. = .~.-countryname+leftright_mid_lag)
lm_mid_ldv_se <- chcse(lm_mid_ldv)

## make table ----
stargazer(lm_mid_fe_base, lm_mid_fe_str, lm_mid_fe,
          lm_mid_ldv_base, lm_mid_ldv_str, lm_mid_ldv,
          se = list(lm_mid_fe_base_se, lm_mid_fe_str_se, lm_mid_fe_se,
                    lm_mid_ldv_base_se, lm_mid_ldv_str_se, lm_mid_ldv_se),
          column.labels = c(rep("FE",3), rep("LDV",3)), 
          omit = c("countryname", "eyear"), omit.labels = c("Country FE", "Year FE"), 
          df = F, digits = 2, align = T, style = "ajps", notes.label = "Notes: ", notes.align ="l", 
          star.cutoffs = c(0.05, 0.01, 0.001), notes = "Heteroscedasticity-consistent standard errors clustered by country in brackets.") 


# Figure 5 ----------------------------------------------------------------

## estimate CMEs ----
me_mid_fe <- margins(lm_mid_fe, vcov=vcovHC(lm_mid_fe, type="HC0", cluster = "countryname"),
                     data = droplevels(lm_mid_fe$model), variables = "leftright_med_lag", 
                     at = list(leftright_den_rlog_lag = prediction::seq_range(lm_mid_fe$model[["leftright_den_rlog_lag"]], 10)))
me_mid_ldv <- margins(lm_mid_ldv, vcov=vcovHC(lm_mid_ldv, type="HC0", cluster = "countryname"),
                      data = droplevels(lm_mid_ldv$model), variables = "leftright_med_lag", 
                      at = list(leftright_den_rlog_lag = prediction::seq_range(lm_mid_ldv$model[["leftright_den_rlog_lag"]], 10)))

## Panel A: FE Estimate ----
ggplot(me_mid_fe) +
  geom_histogram(data = lm_mid_fe$model, aes(leftright_den_rlog_lag, ..count../sum(..count..)*10), 
                 bins = 50,  color = "grey", fill = "white", alpha = 0.7) + 
  geom_line(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag)) +
  geom_line(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag + 1.96*sqrt(Var_dydx_leftright_med_lag)), lty = "dashed") +
  geom_line(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag - 1.96*sqrt(Var_dydx_leftright_med_lag)), lty = "dashed") +  
  geom_point(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag), shape = 18, size = 2) +
  labs(title = "CME of median voter position on main party average position", subtitle = "Conditional on median voter certainty, FE estimates",
       x = "Lagged Median Voter Position Certainty", y = "CME of Lagged Median Voter Position") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10))

## Panel B: LDV Estimate ----
ggplot(me_mid_ldv) +
  geom_histogram(data = lm_mid_ldv$model, aes(leftright_den_rlog_lag, ..count../sum(..count..)*10), 
                 bins = 50,  color = "grey", fill = "white", alpha = 0.7) + 
  geom_line(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag)) +
  geom_line(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag + 1.96*sqrt(Var_dydx_leftright_med_lag)), lty = "dashed") +
  geom_point(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag), shape = 18, size = 2) +
  geom_line(aes(x = leftright_den_rlog_lag, y = dydx_leftright_med_lag - 1.96*sqrt(Var_dydx_leftright_med_lag)), lty = "dashed") +  
  labs(title = "CME of median voter position on main party average position", subtitle = "Conditional on median voter certainty, LDV estimates",
       x = "Lagged Median Voter Position Certainty", y = "CME of Lagged Median Voter Position") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10))

# Table 2 -----------------------------------------------------------------

## define helpers ----
fit_models <- function(.data, .type = c("FE", "LDV")) {
  # make formulas (note diff. for ldvs)
  fmls <- mapply(FUN = function(lhs,rhs) paste0(lhs, "~", rhs), c("lhs_ml","lhs_mr"), paste0(c("rhs_ml","rhs_mr"),"*leftright_den_rlog_lag") )
  if (.type=="FE") fmls <- lapply(fmls, function(fml) as.formula(paste0(fml, "+countryname+eyear")))
  if (.type=="LDV") fmls <- mapply(function(fml, ldv) as.formula(paste0(fml, "+", ldv, "+eyear")), fmls, c("lhs_ml_lag","lhs_mr_lag") )
  # fit models (stepwise)
  m1s <- lapply(fmls, function(fml) lm(data = .data, formula = fml))
  m2s <- lapply(m1s, function(m) update(m, formula. = .~.+gov_party_y3+ud_ipol_y3+openc_y3+realgdpgr_y3+unemp_y3))
  m3s <- lapply(m2s, function(m) update(m, formula. = .~.+compl_leftright+compr_leftright+vturn_lag+leftright_pl_lag+effpar_ele_lag))
  # correct ses 
  correct_ses <- function(m) coeftest(m, vcov. = vcovHC(m, type = "HC0", cluster = "countryname"))[,2]
  cse1s <- lapply(m1s, correct_ses)
  cse2s <- lapply(m2s, correct_ses)
  cse3s <- lapply(m3s, correct_ses)
  # return fits
  ms <- list(m1s,m2s,m3s)
  ses <- list(cse1s,cse2s,cse3s)
  list(
    "ml" = list("lm" = lapply(ms, function(ls) ls[[1]]), "se" = lapply(ses, function(ls) ls[[1]])),
    "mr" = list("lm" = lapply(ms, function(ls) ls[[2]]), "se" = lapply(ses, function(ls) ls[[2]]))
  )
}
export_regtable <- function(ls) { # ls = output from `fit_models`
  stargazer(
    # models 
    c(ls[["ml"]][["lm"]], ls[["mr"]][["lm"]]), 
    # ses
    se = c(ls[["ml"]][["se"]], ls[["mr"]][["se"]]),
    # lab(y)
    dep.var.labels = c("Main Left", "Main Right"), dep.var.labels.include = TRUE,
    # drop lab(FE)
    omit = c("countryname", "eyear"), omit.labels = c("Country FE", "Year FE"), 
    # other aes
    df = F, digits = 2, align = T, style = "ajps", notes.label = "Notes: ", notes.align ="l",
    star.cutoffs = c(0.05, 0.01, 0.001), notes = "Heteroscedasticity-consistent standard errors clustered by country in brackets.") 
}

## fit models ----
lms.fe <- fit_models(.data = data_core2, .type = "FE") # cols 1-6
# lms.ldv <- fit_models(.data = data_core2, .type = "LDV") # not used

## make table ----
export_regtable(ls = lms.fe) # use FE only; manually combine left and right rows in Latex editor

# Figure 6 ----------------------------------------------------------------

## define helpers ----
### get CMEs from model objs
get_cmes <- function(ls) { # ls = output from `fit_models`
  # inner fn
  get_margs <- function(M, XX, XX2 = "leftright_den_rlog_lag") {
    at_XX2 = list(prediction::seq_range(M$model[[XX2]], 10)) |> setNames(XX2) 
    margins(M, vcov=vcovHC(M, type="HC0", cluster = "countryname"), data = droplevels(M$model), variables = XX, at = at_XX2)
  }
  # calc cmes
  m1 = ls[["ml"]][["lm"]][[3]]
  m2 = ls[["mr"]][["lm"]][[3]]
  cme_l = get_margs(M = m1, XX = "rhs_ml")
  cme_r = get_margs(M = m2, XX = "rhs_mr")
  list("ml" = cme_l, "mr" = cme_r)
}
### plot CMEs with CIs
plot_cmes <- function(df_cme, ls_m) { # df_cme = output from `get_cmes`; ls_m = output from `fit_models`
  df <- df_cme %>% 
    map(~{summary(.x) %>% as.data.frame}) %>% 
    `names<-`(c("Main Left","Main Right")) %>% 
    bind_rows(.id = "party")
  df_x <- ls_m %>% 
    map(~{.x[["lm"]][[3]]$model}) %>% 
    map(~{.x %>% dplyr::select(leftright_den_rlog_lag)}) %>% 
    `names<-`(c("Main Left","Main Right")) %>% 
    bind_rows(.id = "party")
  ggplot(df, aes(x=leftright_den_rlog_lag, y=AME, color=party, fill=party, shape=party)) +
    geom_hline(yintercept = 0, lty = "dashed", size = 0.5) +
    geom_ribbon(aes(ymin=lower,ymax=upper), alpha = 0.5, color = NA) +
    geom_line() + geom_point() +
    geom_rug(data = df_x, aes(x=leftright_den_rlog_lag), inherit.aes = FALSE, sides = "b", alpha = 0.5, size = 0.3) +
    facet_grid(cols = vars(party)) +
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    scale_color_grey(start = 0.6, end = 0.2) + scale_fill_grey(start = 0.6, end = 0.2) +
    theme_minimal(base_size = 13) +
    theme(legend.position = "none", strip.text = element_text(size = 14)) +
    labs(
      x = "Lagged Median Voter Position Certainty", y = "CME of Party-Voter Position Difference", 
      title = "CME of party-voter position difference on each party's position change", subtitle = "Conditional on median voter certainty, FE estimates")
}

## estimate CMEs ----
cmes.fe <- get_cmes(lms.fe) # use FE only

## Panels Left & Right ----
plot_cmes(cmes.fe, lms.fe)

# Table 3 -----------------------------------------------------------------

## fit models ----
### cols 1-3: FE models
#### col 1
lm_dif_fe_base <- lm(data = data_core, leftright_dif~leftright_den_rlog_lag+countryname+eyear)
lm_dif_fe_base_se <- chcse(lm_dif_fe_base)
#### col 2
lm_dif_fe_str <- update(lm_dif_fe_base, formula. = .~.+gov_party_y3+ud_ipol_y3+openc_y3+realgdpgr_y3+unemp_y3)
lm_dif_fe_str_se <- chcse(lm_dif_fe_str)
#### col 3
lm_dif_fe <- update(lm_dif_fe_base, 
                    formula. = .~.+
                      gov_party_y3+ud_ipol_y3+openc_y3+realgdpgr_y3+unemp_y3+
                      compl_leftright+compr_leftright+vturn_lag+leftright_pl_lag+effpar_ele_lag)
lm_dif_fe_se <- chcse(lm_dif_fe)
### cols 4-6: LDV models
### col 4
lm_dif_ldv_base <- lm(data = data_core, leftright_dif~leftright_den_rlog_lag+leftright_dif_lag+eyear)
lm_dif_ldv_base_se <- chcse(lm_dif_ldv_base)
### col 5
lm_dif_ldv_str <- update(lm_dif_fe_str, formula. = .~.-countryname+leftright_dif_lag)
lm_dif_ldv_str_se <- chcse(lm_dif_ldv_str)
# col 6
lm_dif_ldv <- update(lm_dif_fe, formula. = .~.-countryname+leftright_dif_lag)
lm_dif_ldv_se <- chcse(lm_dif_ldv)

## make table ----
stargazer(lm_dif_fe_base, lm_dif_fe_str, lm_dif_fe,
          lm_dif_ldv_base, lm_dif_ldv_str, lm_dif_ldv,
          se = list(lm_dif_fe_base_se, lm_dif_fe_str_se, lm_dif_fe_se,
                    lm_dif_ldv_base_se, lm_dif_ldv_str_se, lm_dif_ldv_se),
          column.labels = c(rep("FE",3), rep("LDV",3)), 
          omit = c("countryname", "eyear"), omit.labels = c("Country FE", "Year FE"), 
          df = F, digits = 2, align = T, style = "ajps", notes.label = "Notes: ", notes.align ="l", 
          star.cutoffs = c(0.05, 0.01, 0.001), notes = "Heteroscedasticity-consistent standard errors clustered by country in brackets.")


# save outputs (model objs) -----------------------------------------------

save(list=(ls(pattern = "lm_(mid|dif)")), file = "RData_result_main.RData")

# reset graphx ------------------------------------------------------------

theme_set(theme_base)

# clean env. --------------------------------------------------------------

rm(list = ls())

# END SCRIPT